Analysis
Perform differential gene expression testing
testTreat <- unique(sce$Treatment)
de.res <- mclapply("DMSO", mc.cores = ncores,
function(z) {
message(paste0("Fitting model for ", z))
compDrug <- testTreat[testTreat != z] # instead of subtypes we use other treatments here
mclapply(compDrug, mc.cores = ncores,
function(x) {
message(paste0("Drug ", x))
summed <- sce#[hvgenes, ]
colData(summed)$Treatment <- factor(colData(summed)$Treatment, levels = c(z, compDrug))
y <- DGEList(counts(summed), samples = colData(summed))
design <- model.matrix(~0 + factor(patientID) + Treatment, y$samples) # set coefficient (name or matrix column) - fig. 7 f1000
keep <- filterByExpr(y, design = design)
y <- y[keep,]
summary(keep)
message(paste0("Keeping ", sum(keep), " genes"))
y <- calcNormFactors(y)
y <- estimateDisp(y, design) # gene-wise dispersion (variance on top of poisson distribution)
plotBCV(y)
fit <- glmQLFit(y, design, robust = TRUE) # QL way of fitting the model
plotQLDisp(fit)
res <- glmQLFTest(fit, coef = paste0("Treatment", x))
summary(decideTests(res))
resTab <- res$table
resTab$gene <- rownames(resTab)
resTab <- resTab %>% mutate(p.adj = p.adjust(PValue, method = "BH")) #%>%
#mutate(p.adj = ifelse(p.adj < 1e-12, 1e-12, p.adj))
resTab$Treatment <- x
resTab$Control <- z
resTab %>% arrange(desc(logFC))
}) %>% bind_rows()
}) %>% bind_rows()
## Fitting model for DMSO
de.res <- de.res %>%
left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID"))
de.res %>% filter(Treatment == "Selinexor", Symbol %in% c("BCL2", "TP53"))
## logFC logCPM F PValue gene p.adj
## 1 0.99175500 6.148050 148.6069142 8.633760e-22 ENSG00000141510 6.892299e-19
## 2 0.07380089 6.638112 0.2161808 6.429888e-01 ENSG00000171791 8.493301e-01
## Treatment Control Symbol
## 1 Selinexor DMSO TP53
## 2 Selinexor DMSO BCL2
#write.csv(de.res, opt$deres)
Show the number of DEG
sigTab <- de.res %>%
group_by(Treatment) %>%
filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) %>%
mutate(dir = ifelse(logFC > 0, "Up", "Down")) %>%
dplyr::count(dir) %>%
dplyr::rename(Direction = dir) %>%
ungroup()
drugOrder <- sigTab %>%
group_by(Treatment) %>%
summarise(allDEG = sum(n)) %>%
arrange(desc(allDEG)) %>%
pull(Treatment)
sigTab <- sigTab %>%
mutate(Treatment = factor(Treatment, levels = drugOrder),
n = ifelse(Direction == "Down", n * -1, n))
p <- sigTab %>%
ggplot(aes(x = Treatment, y = n, fill = Direction)) +
geom_bar(stat = "identity", colour = "black") +
geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Up", ], size = 5, nudge_y = 200) +
geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Down", ], size = 5, nudge_y = -200) +
xlab("") + ylab("Number of Genes") +
ggtitle("Significantly Differentially Expressed Genes") +
lgd + scale_fill_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
plot.title = element_text(size = 21, vjust = 1),
axis.text = element_text(size = 17.5),
text = element_text(size = 22.5),
panel.border = element_rect(linewidth = 1))
p
ggsave(plot = p, filename = paste0(opt$plot, "n_deg.png"),
height = 17.5, width = 22, units = "cm")
Construct a network of differentially expressed genes
# Assign nodes (treatments) and edges (5% FDR DEGs)
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25)
allTreat <- unique(plotMat$Treatment)
target <- data.frame(Treatment = allTreat,
target = c("Autophagyi", "IAPi", "HDACi", "mTORi", "ITKi", "TLR8a", "MDM2i",
"IMID", "JAKi", "XPOi", "BCL2i"))
opt$sub <- FALSE
if(opt$sub == TRUE) {
plotMat <- de.res %>% filter(p.adj < 0.1, !Treatment %in% c("Pomalidomide", "Dacinostat")
) %>%
filter(logFC > 0.25 | logFC < -0.25) #%>%
#filter(logFC < -0.25)
} else {
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) #%>%
#filter(logFC > 0.25)
}
allTreat <- unique(plotMat$Treatment)
nodes <- data.frame(Treatment = allTreat)
edges <- plotMat %>% select(Treatment, Symbol) %>%
mutate(Treatment = as.character(Treatment))
# Create network object
set.seed(100)
net <- network(edges, directed = TRUE)
netTab <- ggnet2(net)$data |>
mutate(Type = ifelse(label %in% allTreat, "Treatment", "Symbol"),
Treatment = ifelse(Type == "Treatment", label, NA)) |>
select(-alpha,-color,-shape,-size) |>
mutate(x = as.numeric(x), y = as.numeric(y))
# Define edges based on direction and padj of DEGs
from <- edges %>%
left_join(netTab, by = c("Treatment" = "label")) %>% dplyr::select(Treatment, x, y) %>%
dplyr::rename(x_start = x, y_start = y) |> distinct()
to <- edges %>%
left_join(netTab, by = c("Symbol" = "label")) %>% dplyr::select(Symbol,x,y) %>%
dplyr::rename(x_end = x, y_end = y) |> distinct()
edges <- edges %>% left_join(from) %>% left_join(to) %>%
left_join(plotMat %>% select(Symbol,Treatment,logFC,p.adj)) %>%
distinct() %>%
mutate(Direction = ifelse(logFC < 0,"Down","Up"),
alpha = -log10(p.adj)
#alpha = rescale(-log10(p.adj#squish(p.adj, c(1.812407e-209,1))))
)
# Create auxiliary table for annotation
n_DEG <- plotMat %>% #filter(!Treatment %in% c("Pomalidomide", "BafilomycinA1", "Dacinostat")) %>%
group_by(Treatment) %>%
filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) %>%
dplyr::count() %>% ungroup()
netTabAnno <- netTab |>
left_join(target) |>
left_join(n_DEG) |>
mutate(Treatment = factor(Treatment, levels = allTreat))
p <- ggplot() +
geom_curve(data = edges,
aes(x = x_start, xend = x_end, y = y_start, yend = y_end, col = Direction, alpha = p.adj),
linewidth = 0.1, curvature = 0.05) +
geom_point(data = netTabAnno[netTabAnno$Type == "Symbol",],
aes(x = x, y = y),
size = 0.5, shape = 21, color = "grey60", fill = "grey95", stroke = 0.75, alpha = 0.6) +
geom_point(data = netTabAnno[netTabAnno$Type == "Treatment",],
aes(x = x, y = y, fill = Treatment, size = n), shape = 21, stroke = 0.75) +
# scale_fill_manual(values = col_TX) +
scale_size_continuous(name = "Number of DEGs",
breaks = seq(0, max(n_DEG$n), by = 1000), limits = c(0, max(n_DEG$n)), range = c(0.5,7.5)) +
scale_color_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
scale_alpha_continuous() +
geom_label_repel(data = netTabAnno |> filter(Type == "Treatment"),
aes(x = x, y = y, label = Treatment, fill = Treatment, point.size = rescale(n, to = c(0.5,7.5))),
max.overlaps = Inf, vjust = "bottom", nudge_y = -0.001, fontface = "bold", alpha = 0.7, #direction = "y",
show.legend = FALSE, size = 5) +
guides(color = guide_legend(ncol = 2, override.aes = list(linewidth=2.5)),
size = guide_legend(ncol = 2),
alpha = guide_legend(override.aes = list(linewidth=2.5)),
fill = guide_legend(override.aes = list(size=3.5))) +
xlab("") + ylab("") + ggtitle("Network of Differentially Expressed Genes") +
lgd +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
panel.border = element_rect(linewidth = 1, colour = NA, fill = "transparent")) +
scale_fill_manual(values = c("Venetoclax" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000",
"Ibrutinib" = "#283593", "Everolimus" = "#43A047", "Selinexor" = "#FFD45F",
"BafilomycinA1" = "#F8BBD0", "Motolimod" = "#AB47BC", "Dacinostat" = "#F44336",
"Dacinostat" = "#64B5F6", "Pomalidomide" = "#80DEEA", "Ruxolitinib" = "#4DB6AC"))
p
if(opt$sub == TRUE) {
ggsave(plot = p, filename = paste0(opt$plot, "deg_subnetwork.png"),
height = 20, width = 25, units = "cm")
} else {
ggsave(plot = p, filename = paste0(opt$plot, "deg_network.png"),
height = 20, width = 25, units = "cm")
}
Volcano plots
pList <- lapply(c("Birinapant", "Selinexor", "Nutlin-3a"), function(x) {
if(x %in% c("Motolimod", "Everolimus", "Ibrutinib", "BaflimoycinA1")) {
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
} else if(x %in% c("Birinapant")) {
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF")) %>% select(gene) %>% unlist() %>% as.vector()
} else if(x %in% c("Nutlin-3a", "Selinexor")) {
geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()
geneSelec <- c("TP53", geneSelec)
} else if(x %in% c("Everolimus")) {
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
}
volTab <- de.res %>% filter(Treatment == x)
#volTab <- de.res.2020 %>% filter(Treatment == x, celltype == "T-PLL") %>% dplyr::rename(Symbol = gene)
## Define limits
minLim <- -1 * max(abs(de.res$logFC))
maxLim <- max(abs(de.res$logFC))
thresh <- 0.25
volTab$diffCol <- NA
volTab[volTab$logFC > 0 & volTab$p.adj < 0.1 & volTab$logFC > thresh, "diffCol"] <- "sigUp"
volTab[volTab$logFC < 0 & volTab$p.adj < 0.1 & volTab$logFC < -thresh, "diffCol"] <- "sigDown"
volTab[volTab$p.adj >= 0.1, "diffCol"] <- "ns"
volTab[volTab$logFC < thresh & volTab$logFC > -(thresh), "diffCol"] <- "ns"
p <- volTab %>%
mutate(logFC = ifelse(logFC > maxLim, maxLim, logFC)) %>%
mutate(logFC = ifelse(logFC < minLim, minLim, logFC)) %>%
#filter(gene %in% geneSelec) %>%
filter(Treatment == x) %>%
ggplot(aes(x = logFC, y = -log10(p.adj), fill = diffCol)) +
geom_point(shape = 21, size = 2) + xlab("Log2FC") + ylab("-Log10(p.adj)") + ggtitle(paste0(x)) +
#xlim(minLim, maxLim) +
geom_hline(yintercept = -log10(0.1), linetype = "dashed") +
geom_vline(xintercept = -thresh, linetype = "dashed") + geom_vline(xintercept = thresh, linetype = "dashed") +
geom_label_repel(data = filter(volTab, Symbol %in% c(geneSelec), Treatment %in% x, p.adj < 0.1,
logFC > thresh | logFC < -(thresh)), aes(label = Symbol), fill = "white",
nudge_y = 0.25, size = 4, max.overlaps = getOption("ggrepel.max.overlaps", default = 30),
segment.linetype = 2, force = 20, alpha = 0.8) +
scale_fill_manual(values = c("sigUp" = "#F44336", "sigDown" = "#1976D2", "ns" = "grey")) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 17.5),
axis.text = element_text(size = 15),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.background = element_rect(fill='transparent'),
strip.background = element_blank())
ggsave(plot = p, filename = paste0(opt$plot, "vol_", x, ".png"),
height = 15, width = 15, unit = "cm")
p
})
## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
pList
## [[1]]
##
## [[2]]
## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## [[3]]
## Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
geneSub <- lapply(unique(de.res$Treatment), function(x) {
de.res %>%
filter(Treatment == x) %>%
filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) %>%
pull(Symbol)
}) %>% unlist() %>% unique()
geneSub <- intersect(geneSelec, geneSub)
plotMat.wide <- de.res %>%
mutate(logFC = ifelse(p.adj >= 0.1, 0, logFC)) %>%
dplyr::select(Symbol, Treatment, logFC) %>%
filter(Symbol %in% geneSub) %>%
dplyr::rename(log2FC = logFC) %>%
pivot_wider(names_from = "Treatment", values_from = "log2FC") %>%
column_to_rownames("Symbol")
geneOrder <- rownames(plotMat.wide)[hclust(dist(plotMat.wide), method = "ward.D2")$order]
Dot plot
p <- de.res %>% filter(#Treatment %in% c("Birinapant", "Motolimod"),
Symbol %in% geneSub) %>%
filter(p.adj < 0.1) %>%
mutate(pSign = -log10(p.adj) * sign(logFC),
p.adj = ifelse(-log10(p.adj) > 5, 5, -log10(p.adj)),
Symbol = factor(Symbol, levels = geneOrder)) %>%
dplyr::rename(log2FC = logFC) %>%
ggplot(aes(x = Symbol, y = Treatment)) +
geom_point(aes(size = p.adj, fill = log2FC), shape = 21) +
xlab("") + ylab("") + ggtitle("Effect on TNF-NFkB and TLR Target Genes") +
lgd +
scale_fill_gradient2(high = "#F44336", low = "#1976D2", name = "log2FC",
limits=c(-2, 2), oob=squish) +
scale_size_continuous(name = "-Log10(p.adj)",
breaks = seq(0, 5, by = 2.5),
limits = c(0, 5), #range = c(0.5,7.5)
) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
text = element_text(size = 12.5),
panel.border = element_rect(linewidth = 1))
p
ggsave(plot = p, filename = paste0(opt$plot, "biri_moto_nfkb.png"),
height = 9, width = 50, units = "cm")
GSEA
go.df.merge <- mclapply(c("Birinapant", "Motolimod", "Selinexor", "BafilomycinA1", "Dacinostat", "Pomalidomide", "Ibrutinib",
"Everolimus", "Venetoclax", "Ruxolitinib", "Nutlin-3a"), mc.cores = ncores,
function(compTreat) {
message(paste("Running gprofiler for", compTreat, sep = " "))
sigUp <- de.res %>% filter(logFC > 0.25, Treatment == compTreat,
#!gene %in% c(drugOther),
p.adj <= 0.1) %>% arrange(desc(logFC))
#sigUp <- de.res %>% filter(celltype == compType, logFC > 0.25, Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
sigUp
if(nrow(sigUp) > 1) {
gostresUp <- gost(query = sigUp$gene,
organism = "hsapiens", ordered_query = TRUE,
multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.1, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
gostresUp$result
go.df.up <- gostresUp$result %>% arrange(p_value)
go.df.up$rank <- 1:nrow(go.df.up)
go.df.up$type <- "Up"
} else {
go.df.up <- data.frame(p_value = NA, rank = NA, type = "Up", source = NA, term_name = NA,
intersection_size = 10)
}
sigDown <- de.res %>% filter(logFC < -(0.25), Treatment == compTreat,
#!gene %in% c(drugOther),
p.adj <= 0.1) %>% arrange(logFC)
#sigDown <- de.res %>% filter(celltype == compType, logFC < -(0.25), Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
if(nrow(sigDown) > 1) {
gostresDown <- gost(query = sigDown$gene,
organism = "hsapiens", ordered_query = TRUE,
multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.1, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
gostresDown$result
go.df.down <- gostresDown$result %>% arrange(p_value)
go.df.down$rank <- 1:nrow(go.df.down) * -1
go.df.down$type <- "Down"
} else {
go.df.down <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA,
intersection_size = 10)
}
go.df.merge <- rbind(go.df.up, go.df.down)
go.df.merge$Treatment <- compTreat
go.df.merge <- go.df.merge %>% filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(p.adj = p.adjust(p_value, method = "BH"))
if(nrow(go.df.merge) < 1) {
go.df.merge <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA,
intersection_size = 0)
} else {
go.df.merge
}
}) %>% bind_rows()
Show number of enriched gene-sets
#### Show number of enriched gene sets
sigTab <- go.df.merge %>%
dplyr::count(Treatment, type) %>%
dplyr::rename(Direction = type) %>%
mutate(n = ifelse(Treatment %in% c("Ruxolitinib", "Venetoclax"), 0, n)) %>%
mutate(Treatment = factor(Treatment, levels = drugOrder),
n = ifelse(Direction == "Down", n * -1, n)) %>%
filter(!is.na(Treatment))
naTab <- data.frame(Treatment = c("Venetoclax", "Venetoclax", "Selinexor"),
Direction = c("Down", "Up", "Down"),
n = c(0, 0, 0))
sigTab <- rbind(sigTab, naTab)
p <- sigTab %>%
ggplot(aes(x = Treatment, y = n, fill = Direction)) +
geom_bar(stat = "identity", colour = "black") +
geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Up", ], size = 5, nudge_y = 200) +
geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Down", ], size = 5, nudge_y = -200) +
xlab("") + ylab("Number of Gene Sets") +
ggtitle("Significantly Enriched Gene Sets") +
lgd + scale_fill_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
plot.title = element_text(size = 21, vjust = 1),
axis.text = element_text(size = 17.5),
text = element_text(size = 22.5),
panel.border = element_rect(linewidth = 1))
p
ggsave(plot = p, filename = paste0(opt$plot, "n_gsea.png"),
height = 17.5, width = 22, units = "cm")
Make dot plot
## Exclude all terms that are up and downregulated by a single-drug
dblTerm <- go.df.merge %>%
filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
#filter(term_name %in% unique(termTab$term_name)) %>%
dplyr::count(term_name, Treatment) %>%
filter(n > 1) %>%
pull(term_name) %>% unique()
## Keep the top20 of each compound
nTop <- 5
termTab <- lapply(c("Selinexor", "Nutlin-3a", "Venetoclax", "Birinapant", "Ibrutinib", "Motolimod",
"Dacinostat", "Pomalidomide", "Everolimus"), function(y) {
termTab <- lapply(c("Up", "Down"), function(z) {
message(y)
message(z)
gseaTab <- go.df.merge %>% filter(type == z, Treatment == y) %>%
filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
filter(!term_name %in% dblTerm)
gseaTab <- gseaTab[1:nTop, ] # get only top enrichments
#termOrder <- gseaTab %>% arrange(desc(p_value)) %>% pull(term_name)
gseaTab <- gseaTab %>%
mutate(p.adj = p.adjust(p_value, method = "BH"))
}) %>% bind_rows()
}) %>% bind_rows()
gseaTab <- go.df.merge %>%
filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
filter(term_name %in% unique(termTab$term_name)) %>%
filter(!term_name %in% dblTerm)
plotMat <- gseaTab %>%
mutate(p.adj = -log10(p.adj), p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
select(term_name, Treatment, p.adj) %>% unique() %>%
pivot_wider(names_from = "Treatment", values_from = "p.adj") %>%
column_to_rownames("term_name")
plotMat[is.na(plotMat)] <- 0
geneOrder <- rownames(plotMat)[hclust(dist(plotMat), method = "ward.D2")$order]
p <- gseaTab %>%
mutate(term_name = factor(term_name, levels = geneOrder),
intersection_size = ifelse(intersection_size > 200, 200, intersection_size)) %>%
mutate(p.adj = -log10(p.adj), p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
ggplot(aes(x = term_name, y = Treatment)) +
geom_point(aes(size = intersection_size, fill = p.adj), shape = 21) +
xlab("") + ylab("") + ggtitle("Top Enriched Gene Sets") +
lgd +
theme(axis.text = element_text(size = 10),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 8.5),
panel.background = element_rect(linewidth = 1, fill = "transparent")) +
scale_size_continuous(name = "Number of \nOverlapping Genes",
breaks = seq(0, max(gseaTab$intersection_size), by = 50),
limits = c(1, 200), range = c(0.5,7.5)) +
scale_fill_gradient2(low = "#1976D2", high = "#F44336", limits=c(-4, 4), oob=squish,
name = "-Log10(p.adj) \nwith Direction")
p
ggsave(plot = p, filename = paste0(opt$plot, "top_genesets_dotplot.png"),
height = 19.5, width = 50, units = "cm")
Output session info
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel grid stats4 stats graphics grDevices datasets
## [8] utils methods base
##
## other attached packages:
## [1] sna_2.7-2 statnet.common_4.9.0
## [3] scales_1.3.0 numbat_1.4.0
## [5] Matrix_1.6-1.1 airr_1.5.0
## [7] ggsci_3.2.0 BiocParallel_1.36.0
## [9] ggnet_0.1.0 network_1.18.2
## [11] tftargets_1.3 ggrepel_0.9.5
## [13] readxl_1.4.3 edgeR_4.0.16
## [15] limma_3.58.1 circlize_0.4.16
## [17] ComplexHeatmap_2.18.0 gridExtra_2.3
## [19] gprofiler2_0.2.3 scater_1.30.1
## [21] scran_1.30.2 scuttle_1.12.0
## [23] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
## [25] Biobase_2.62.0 GenomicRanges_1.54.1
## [27] GenomeInfoDb_1.38.8 IRanges_2.36.0
## [29] S4Vectors_0.40.2 BiocGenerics_0.48.1
## [31] MatrixGenerics_1.14.0 matrixStats_1.3.0
## [33] lubridate_1.9.3 forcats_1.0.0
## [35] stringr_1.5.1 dplyr_1.1.4
## [37] purrr_1.0.2 readr_2.1.5
## [39] tidyr_1.3.1 tibble_3.2.1
## [41] ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 ggplotify_0.1.2
## [3] cellranger_1.1.0 polyclip_1.10-7
## [5] lifecycle_1.0.4 doParallel_1.0.17
## [7] lattice_0.21-9 MASS_7.3-60
## [9] magrittr_2.0.3 plotly_4.10.4
## [11] sass_0.4.9 rmarkdown_2.27
## [13] jquerylib_0.1.4 yaml_2.3.9
## [15] metapod_1.10.1 RColorBrewer_1.1-3
## [17] abind_1.4-5 zlibbioc_1.48.2
## [19] quadprog_1.5-8 hahmmr_1.0.0
## [21] ggraph_2.2.1 RCurl_1.98-1.14
## [23] yulab.utils_0.1.5 tweenr_2.0.3
## [25] GenomeInfoDbData_1.2.11 irlba_2.3.5.1
## [27] tidytree_0.4.6 dqrng_0.4.1
## [29] DelayedMatrixStats_1.24.0 codetools_0.2-19
## [31] DelayedArray_0.28.0 ggforce_0.4.2
## [33] tidyselect_1.2.1 shape_1.4.6.1
## [35] aplot_0.2.3 farver_2.1.2
## [37] ScaledMatrix_1.10.0 viridis_0.6.5
## [39] jsonlite_1.8.8 GetoptLong_1.0.5
## [41] BiocNeighbors_1.20.2 tidygraph_1.3.1
## [43] iterators_1.0.14 systemfonts_1.1.0
## [45] foreach_1.5.2 tools_4.3.2
## [47] ragg_1.3.2 treeio_1.29.0.002
## [49] Rcpp_1.0.12 glue_1.7.0
## [51] SparseArray_1.2.4 xfun_0.45
## [53] withr_3.0.0 fastmap_1.2.0
## [55] bluster_1.12.0 fansi_1.0.6
## [57] digest_0.6.36 rsvd_1.0.5
## [59] parallelDist_0.2.6 timechange_0.3.0
## [61] R6_2.5.1 gridGraphics_0.5-1
## [63] textshaping_0.4.0 colorspace_2.1-0
## [65] RhpcBLASctl_0.23-42 utf8_1.2.4
## [67] generics_0.1.3 renv_1.0.7
## [69] data.table_1.15.4 graphlayouts_1.1.1
## [71] httr_1.4.7 htmlwidgets_1.6.4
## [73] S4Arrays_1.2.1 pkgconfig_2.0.3
## [75] gtable_0.3.5 XVector_0.42.0
## [77] htmltools_0.5.8.1 clue_0.3-65
## [79] png_0.1-8 ggfun_0.1.5
## [81] knitr_1.48 rstudioapi_0.16.0
## [83] tzdb_0.4.0 rjson_0.2.21
## [85] coda_0.19-4.1 nlme_3.1-163
## [87] cachem_1.1.0 GlobalOptions_0.1.2
## [89] vipor_0.4.7 pillar_1.9.0
## [91] logger_0.3.0 vctrs_0.6.5
## [93] BiocSingular_1.18.0 beachmat_2.18.1
## [95] cluster_2.1.4 beeswarm_0.4.0
## [97] evaluate_0.24.0 cli_3.6.3
## [99] locfit_1.5-9.10 compiler_4.3.2
## [101] rlang_1.1.4 crayon_1.5.3
## [103] labeling_0.4.3 fs_1.6.4
## [105] ggbeeswarm_0.7.2 stringi_1.8.4
## [107] viridisLite_0.4.2 munsell_0.5.1
## [109] lazyeval_0.2.2 hms_1.1.3
## [111] patchwork_1.2.0 sparseMatrixStats_1.14.0
## [113] statmod_1.5.0 highr_0.11
## [115] igraph_2.0.3 memoise_2.0.1
## [117] scistreer_1.2.0 RcppParallel_5.1.8
## [119] bslib_0.7.0 phangorn_2.11.1
## [121] ggtree_3.10.1 fastmatch_1.1-4
## [123] ape_5.8